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ABSTRACT. Here, through a systematic methodology and the use of high performance computing, we calculate the optimum 
shape for a wave energy converter under the action of incident waves of (i) monochromatic unidirectional, (ii) monochro- 
matic directional, (iii) polychromatic unidirectional and (iv) polychromatic directional (with both directional symmetry and 
asymmetry). As a benchmark for our study, without loss of generality, we consider a submerged planar pressure differential 
wave energy converter, and use Genetic Algorithm to search through a wide range of shapes. A new parametric description 
of absorber shape based on Fourier decomposition of geometrical shapes is introduced, and for each shape hydrodynamic 
coefficients are calculated, optimum power take-off parameters are obtained, and overall efficiency is determined. We show 
that an optimum geometry of the absorber plate can absorb a significantly higher energy (in some cases few times higher) 
when compared to a circular shape of the same area. Specifically, for a unidirectional incident wave, the optimum shape, 
as expected, is obtained to be the most elongated shape. For directional incident waves, a butterfly- shape is the optimum 
geometry whose details depend on not only the amplitude and direction of incident wave components, but also the relative 
phases of those components. For the latter effect, we find an optimally averaged profile through a statistical analysis. 
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Oceanic surface waves carry a much higher energy density, sometimes by two orders of magnitude, than wind 
and solar. Estimates show that ocean wave energy can realistically provide about 10% of world’s electricity need 
[1-6]. This potential, along with proximity to load centers that are typically along coastal areas, as well as good 
predictability and power consistency makes wave energy an appealing solution to our power-hungry world that is 
fed up with archaic pollutant-rich energy production methods [7-10]. Harnessing ocean wave energy, however, faces 
several major challenges that have, so far, barred wave energy industry from taking off. These include complexities 
associated with working in the harsh ocean environment such as corrosive nature of the ocean water, or extreme loads 
and impacts during storms. Even in a typical day, incoming ocean waves arrive almost always as a spectrum, composed 
of waves of different frequencies and directions. For a wave energy device which is required to work at resonance, 
this poses an extra challenge; which is also new to the industry and research: conventionally offshore structures are 
designed to avoid resonance as much as possible, and to stay put by minimizing all motions. A wave energy device on 
the other hand must target resonance and maximize its motion. 

Over the past century, more than one thousand patents have been filed on different ideas for harvesting ocean wave 
power. Over the past few decades, numerous clever theoretical, computational and experimental advancements have 
resulted in a much clearer picture of what can be achieved from ocean waves [cf. e.g. 11-16]. Most of our today’s 
knowledge, however, is for the case of monochromatic unidirectional incident waves, while ocean waves come in 
spectrum of frequencies and from many different directions. What makes the investigation of a frequency-spectrum 
difficult is frequency-dependent coefficients in the governing equation of a wave energy device that prevents closed- 
form solutions. High-performance computational tools are now allowing computational study of a wide range of 
scenarios that happen under the action of broadband frequency and directional waves. 


E-mail address: soes@stanford.edu, reza.alam@berkeley.edu. 
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Nomenclature 

PTO power take-off unit 

A 

wave length [m] 

BEM boundary element method 

d 

absorber plate’s thickness [m] 

CWR 

. capture width ratio 

e 

wave’s angle of incidence [°] 


shape difference factor 

(Op 

peak wave’s angular frequency 

hopr 

absorber plate’s operating depth [m] 

r 

peak enhancement factor 

hbtm 

sea bottom depth [m] 

a p 

Philips parameter 

Z 

vertical (heave) displacement [m] 

H s 

significant wave height 

K 

spring’s coefficient [N/m] 

D(e ) 

spreading function 

C 

damper’s coefficient [Ns/m] 

I'm 

wave’s mean angle of incidence [°] 

a 

wave amplitude [m] 

© 

directional spectrum’s angular span [°] 

m 

absorber plate’s mass [£g] 


asymmetry parameter 

A 

added mass [£g] 

S (ft), 0)wave’s spectral density 

B 

added damping [Ns / m] 

ro 

radius of base circle [m] 

CO 

wave’s angular frequency [rad /s] 

A 

area of a shape [m 2 ] 

1 exc 

wave induced excitation force [A] 

a n 

coefficient of n th order Fourier mode [m] 

<t>F 

phase of excitation force [rad] 


phase of n th order Fourier mode [rad] 


phase of system response [rad] 

N c 

total number of Fourier modes 

□ 

amplitude of □ 

PA 

elongation coefficients 

r e 

radius of equivalent circle 

N e 

number of directional spreading subdivisions 

T 

wave period [s] 

Nco 

number of frequency spreading subdivisions 

kp 

peak wave number [rad /m] 

P, 

normalized time-averaged extracted power 

k 

wave number [rad /m] 

N 

Genetic Algorithm’s number of generations 

Sf(w 

frequency spectrum 

Pabs 

time-averaged extracted power by the optimum shape [W] 

T P 

peak wave period [s] 

P . 

1 CUT 

time-averaged extracted power by the equivalent circle [W] 


To maximize the efficiency of a wave energy converter, a number of strategies are employed. For instance, we know 
today that (i) depending on the wave condition, depth of the water and the bathymetry, a different class of wave energy 
converters may be more suitable (i.e. more efficient and cost effective) [17-22], (ii) A wide range of active control 
strategies have been proposed and their performance have been investigated [23-30], (iii) Nonlinear mechanisms have 
been shown to enhance the bandwidth of high performance of wave energy devices [31, 32], (iv) Flexibility of the 
absorber may positively contribute (in some cases significantly) to the efficiency of wave energy harnessing devices 
[33], and (v) optimized geometry, shape, and power take-off mechanism of a wave energy converter play an important 
role in the overall performance of the device [34-37], Shape optimization is of relevance in all classes of wave energy 
conversion techniques such as oscillating water columns where it is found that immersion depth of 0.45 of water depth 
and chamber diameter of 0.92 of water depth yields the maximum efficiency [38]; and overtopping devices where it is 
found that the optimum length to opening ratio has to he ~ 2.5 — 3 [39]. Most efforts so far, nevertheless, have been 
limited to comparing few specific shapes. As mentioned above, this has been mainly due to the computational costs 
of finer searches among all possibilities. 

For a heaving point absorber wave energy converter under wave conditions of Belgian coastal area of the North 
Sea (shallow water), Vantorre et al. [40] compared performance of a set of geometries including a hemisphere and a 
few conical geometries. They found that a 90° cone with cylindrical extension that pierces the water surface yields 
the largest efficiency among their set. For a point absorber working in wave conditions off the west coast of Ireland, 
Goggins and Finnegan [41] considered variations of a vertical cylinder (e.g. truncated cylinder with a hemispherical 
base, with a 45° linear chamfer, or with a quarter encircle chamfer), and found that a truncated cylinder with a 
hemisphere attached to its base with the total draft to radius ratio of 2.5 is the optimum shape having the largest 
significant heave velocity response. Through experimentations, Hager et al. [42] showed that a buoy with concave 
faces performs better than a buoy with flat or convex faces, and that the capture width ratio can reach up to 94.5%. For 
a singe motion of a wave energy device, McCabe [43] parametrized the device geometry by bi-cubic B-spline surfaces 
and showed that the optimum shape has asymmetry in the direction of wave propagation with distinct, pointed prows 
and sterns. 

Here through a systematic approach, we find optimum shape of a single wave energy converter under different 
wave conditions including monochromatic, polychromatic, unidirectional and directional seas. While the methodology 
proposed is general and applies to any wave energy converter design, we focus our attention on the optimization of 
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Wave Direction 



Sea Bottom 


Figure 1 . Schematic of a submerged planar pressure differential wave energy converter in sea; 
h opr , a , h/,i m , K , C are absorber plate’s operating depth, wave amplitude, sea bottom depth, and 
power take-off’s spring and damper coefficients respectively 


a submerged planar pressure differential wave energy converter called “ Wave Carpet" [44 — 46] for deep waters. The 
original idea of the wave carpet was to mimic the natural phenomenon of surface waves getting damped by muddy 
seabeds, and to place a synthetic wave absorber near the seabed that responds to the action of overpassing waves 
similar to how mud responds. If water is deep, the device may be elevated in the water column so that the device 
stays close enough to the surface, hence retaining required performance. A submerged wave energy converter is more 
survivable under storm conditions (e.g. resulting slamming impacts [47]) as the water column buffers surface forces, 
does not interfere with surface vessels and fishing boats, has no visual pollution, and does not impede atmosphere-sea 
surface interactions such as Oxygen and C02 exchange, which is of major concerns from environmental point of view 
particularly when large-scale deployments (e.g. wave energy farms [48, 49]) are sought. 

In our investigation, without loss of generality, we consider a single deep water version of the aforementioned wave 
energy converter, a rigid planar absorber (figure 1), and restrict our analysis to its heave-only motion. The objective 
is to find the optimum shape that produces maximum power under the action of specific form of incident waves. In 
our optimization, we compare performance of shapes of the same area. This is because the overall cost of ships and 
offshore structures is roughly proportional to their weights. Therefore, for a wave energy device to have a comparable 
price with another of the same type (say Wave Carpet in this study), they must have, roughly speaking, the same weight 
and therefore the same area. Clearly this is not the only factor determining the cost [cf. 50], and it is not difficult to 
find variations of a shape with the same area whose production leads to costs orders of magnitude higher than other 
shapes (e.g. infinitely long and thin plate). To compensate for such effects, we consider that all permissible shapes are 
variations of a fixed base ellipse plus a few Fourier modes, but put a limit on how much elongated the ellipse can be, 
and also on the amplitude of Fourier modes to avoid sharp corners. For each given wave condition, we optimize (i) 
spring stiffness, (ii) damping coefficient of generator, (iii) major and minor axes of the ellipse, and (iv) amplitude and 
phase of each Fourier mode (the last two items iii and iv under the constraint that the total area is the same). 

As optimization algorithm, we tried different strategies including gradient-based methods, direct search methods, 
stochastic methods, and Genetic Algorithm. It turns out that for the problem in hand Genetic Algorithm has the fastest 
convergence rate, and we confirmed that it in fact converges to the true solution via cross-validating its results with 
those obtained via the direct search methods. Genetic Algorithm is a heuristic search method that has gained attention 
due to its efficient capability in finding optimum solutions; particularly in cases where the objective function cannot 
be expressed analytically in terms of its constitutive parameters/variables and where there exists many local extrema 
in the search area [51]. Genetic algorithm is not foreign to offshore industry, and has been used for optimization of 
offshore structures [e.g. 43, 50, 52], and also for optimizing wave energy devices [e.g. 43, 50]. 

Here, we start with problem statement and presentation of the governing equations (§2), followed by details of our 
shape optimization methodology and definition of the objective function (§3). We then present and discuss results 
of optimization for the absorber’s shape under four different wave conditions (§4): (i) monochromatic unidirectional, 
(ii) monochromatic directional, (iii) polychromatic unidirectional and (iv) polychromatic directional waves (with both 
directional symmetry and asymmetry). 
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1. Fundamentals 


We consider a single submerged flat plate with uniform thickness of d in ’-direction at the depth h opr connected 
to a power take-off (PTO) system and under the action of surface gravity waves (see figure 1). The power take-off 
system is composed of a linear spring with the stiffness K, and a generator modeled as a linear damper with a damping 
coefficient C. The idea is that the net pressure difference between the top and the bottom of the plate due to the passage 
of surface gravity waves can move the plate that in turn moves the end of the power take-off system and results in the 
generation of power. Therefore, energy of surface waves is harnessed by the power take-off unit. 

In our analysis, we assume that water is homogeneous, inviscid, incompressible and irrotational such that potential 
flow theory applies. We also neglect the effect of surface tension, and assume that the water depth is constant and 
equal to hb tm . Restricting the motion of the plate to the heave-direction only (i.e. ’-direction in figure 1 ) the governing 
equation for the motion of absorber plate reads 

(1.1) (m+A)z+(C+B)z + Kz = F exc 


in which m is the mass of the absorber plate, A and B are the hydrodynamic added mass, and added damping coef- 
ficients respectively, and F exc is the wave excitation force in the z-direction. Parameters A . B and F exc are functions 
of the absorber plate’s geometry as well as the incident wave frequency (F exc is also linearly proportional to incident 
wave amplitude). In order to calculate A . B and F exc we use an open source Boundary Element Method-based solver, 
NEMOH, developed at Ecole Centrale de Nante by Babarit and Delhommeau [53] 1 . 

If the incident wave is monochromatic with the frequency ft) then F exc = F e 1C0t where F is the amplitude of the 
excitation. The absorber’s response therefore is obtained readily as z = ze ' (< ™ 1 where 


(1.2a) 


(1.2b) 


z = F/ yj [K — (m + A) ft) 2 ] 2 + ft) 2 [C + B] 2 , 


0 - = arctan 


co(C + B) 

K — (m +A) ft) 2 J 


The absorbed power is therefore 
(1.3) 


r T i 

Pabs=^ I C Z 2 dt = ~CC0 2 Z 2 . 


1 

T . 


'o 


For validation of numerical calculation of hydrodynamic coefficients (added mass, added damping, and excita- 
tion force), we compare the Capture Width Ratio (CWR) 2 of a floating hemisphere with hydrodynamic coefficients 
obtained (i) analytically [58] and (ii) numerically through NEMOH. Figure 2 shows the comparison of results from 
analytical expressions (solid lines) and those obtained numerically (markers), for two cases of parameters tuned such 
that resonance occurs at kr = 1 (case A) and kr = 1.5 (case B). 

The goal of the current research is to look for an optimized shape of a wave energy converter under four different 
sea conditions: i. normal incidence of monochromatic unidirectional waves, ii. monochromatic waves with directional 
spreading (figure 3. a), iii. normal incidence of polychromatic unidirectional waves, and iv. polychromatic directional 
spectrum of waves with a symmetric directional distribution (figure 3.b) and with a general (i.e. asymmetric) direc- 
tional distribution (figure 3.c). For the directional spectra (wave climates ii, iv) the basic premise is that the spectral 
density function S(co,6) is the product of a frequency spectrum Sj ( ft ) ) and a spreading function I) (0). For the fre- 
quency spectrum, we choose the Joint North Sea Wave Project (JONSWAP) spectrum in the form (cf. [59]) 


(1.4) 


c a pS 

Sfico) = — ^exp 

1 ft ) 5 


yexp[-(<B-ffl,fi 2 /(2CT 2 (Bp 2 )] , 


where ft) is the wave frequency, co p the peak wave frequency, and y the peak enhancement factor specifying the 
spectral bandwidth typically ranging between 1 <y<9 for which we choose the mean value of y = 3.3 [59], a p = 
H x a> p /[16/o(y)g 2 ] is the Phillips parameter, H s is significant wave height, and the Zeroth-order moment /o(y) which 
varies in the range of 0.2</o(y)<0.5 is usually calculated numerically [60] being equal to 7 q( 3.3) = 0.3 for our applica- 
tion. In this manuscript, without losing generality we focus on the sea state five (rough sea conditions) with H s = 3.25 
meters and T p = 9.7 seconds [61]. For the directional spreading component of the directional spectra (wave climates ii. 


'Alternative tools include ANSYS AQWA ([40, 54]) and WAMIT ([43, 50, 52, 55-57]). 
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A measure of a wave energy converter’s performance is the Capture Width Ratio (CWR) defined as the power extracted by the wave energy 
converter divided by the amount of power in the incident wave of crest length equal to the width of the device. For instance, if a device of width wj 
captures pd Watt under incident waves of power flux p w Watt per meter of the wave crest, then CWR — pdj ( PwWd )• 
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Figure 2. Capture width ratio (CWR) of a floating hemisphere (with relative density one) in heave 
as a function of dimensionless wavenumber of the incident wave kr ( r is the radius of the hemi- 
sphere). Spring coefficient is chosen such that the undamped system is at resonance at respectively 
kr = 1.0 (Case A), and kr = 1.5 (Case B). Likewise damping coefficient of the power take-off is 
chosen to be the same as the radiation damping for kr = 1, 1.5. Solid lines show CWR calculated 
analytically [58], and markers show results obtained from NEMOH. 



Figure 3. Spectral density function (5) distribution of (a) monochromatic waves with directional 
spreading for four angular spans of ±20°, ±40°, ±60° and ±80°, (b) a JONSWAP broadband 
directional spectrum with symmetric angular distribution, fx = 0,s = 2, 0 = 160°,//, = 3.25 \m],T„ = 
9.7 [.v] , (c) a JONSWAP directional spectrum with an asymmetric angular distribution, fi = 0.5, s = 
75, © = 160° , H s =3.25 [m] , T p =9.7 [.s] 


iv), without loss of generality, we assume zero degree mean for the wave direction (i.e. 9 m = 0°) and use the following 
directional spreading function D(9) (c.f. e.g. [62]) 


(1.5a) 


(1.5b) 


Dm = ( !“»•(§{) e§K0/2 

\ 0 I 0§ I >0/2 

f exp (— /x) 6 > 0 

’ [ exp l+n) 6 <0 


where 0 is angular span, and /x is the asymmetry parameter. Physically speaking, a positive fx means a broader 
angular distribution of energy for 0 > 0 m = 0 (see e.g. figure 3.c). Here, for a symmetric directional spreading 
we choose jx = 0, s = 2, © = 160°, and for asymmetric directional spreading we choose jx = 0.5, © = 160°, and 
s = 75 to limit the spectrum in a finite angular span. These cases are shown in figure 3. For a monochromatic 
directional spectrum, amplitude of individual waves are aj = 2S(6j)d9 , and for a broadband directional spectrum 
aj = y/2 S((Oj, 0j ) d co d 0 . 
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2. Shape Optimization Methodology 


The objective of optimization is to maximize the normalized time-averaged power generation P r defined by 


( 2 . 1 ) 


Pr = 


Pabs 

Pcirc 


where P a \ K is the time-averaged extracted power by a wave energy converter, and P nrr is the time-averaged extracted 
power by a circular-shaped device of the same area (equivalent circle). Note that both Pubs, Parc are obtained for 
optimum power take-off parameters C. K . 

The power extracted for the case of unidirectional monochromatic incident waves (case i) is very simple and given 
by (1.3). But care must be taken in dealing with incident waves with directional spreading. The reason is that when 
we have directional spreading two waves of the same frequency may superimpose at the location of the device, but 
the extracted energy cannot be linearly superimposed as the contribution of phase-differences become of leading order 
importance. 3 Specifically, power absorption under monochromatic directional incident waves (wave climate ii ) is 
given by (see Appendix I) 


( 2 . 2 ) 


j Ng Ng 

Pabs = ~C(0 2 £ £ ZiZj COS - <j> ZJ ) , 
z i= i y= i 


and power absorption under an incident directional spectrum (wave climate tv) is given by 


(2.3) 


2 N (0 ^ Ng Ng 

Pabs = X? E a l E E W COS ^z,i - 

Z k= 1 i=l 7=1 


where z is the system response amplitude, (j) z is the response phase, and Nq,Nco are the total number of directional and 
frequency spreading subdivisions (discretizations). 

For shape parameterization, we consider a horizontal circle of radius ip (the base circle), and add n = N c Fourier 
components with amplitude a„ and phase <j>„ according to 

N c 

(2.4) r(0) = ro + ^ a„ sin [n 9 + 

n= 1 


This representation enables us to construct any arbitrary shape if N c is large enough. Note that with the introduction 
of phases the obtained shapes can be symmetric or asymmetric (cf. exempli gratia [63, 64]). For the problem in hand, 
in some cases, the shape has to be elongated in the x or y-direction, therefore to make the optimization more efficient 
we define elongation coefficients p,q according to 


(2.5) 


x = prcos(d), y = qrsin(Q). 


These elongation coefficients p,q will be optimized along with amplitudes a n and phases to obtain the optimum 
shape. We also impose the constraint that the total area of all shapes must be constant. Mathematically this can be 
expressed, using Green’s theorem [65] (for derivation see Appendix III), as 


(2.6) 


pq 



= Constant. 


For optimization, we use Genetic Algorithm ( GA ) which is an evolutionary optimization scheme that mimics the 
process of natural selection (see e.g. [66]). The GA solver of MATLAB optimization toolbox, Optimtool, is used 
as the optimization module in our work. An optimization through Genetic Algorithm starts with an initial random 
set of points (i.e. candidate solutions) in the phase space called initial population. Fittest members of the initial 
population, i.e. data points that better fulfill the optimization objective, are given preferences to breed and form the 
next generation. To form the next generation, each pair of survived data points, i.e. parents, give birth, through a 
process called crossover to one or more children that inherit properties from both parents, though not being identical 
to either. To make sure that the artificial evolution does not converge to a local extremum, mutation is induced here and 
there through randomly modifying a percentage of newborn children. There are many other details that have shown 


o 

This is similar in nature to the paradoxical classic problem of two same-frequency waves of amplitudes a superposition. Let’s assume we have 
two identical waves of amplitude a, each carry energy of ca 2 with c being a constant. Therefore the two have 2ca 2 energy. If we superpose these 
two waves, there is a new wave of amplitude 2 a that must have, according to the wave theory, the energy of 4 ca 2 which is double the initial energy 
of the two initial waves. Alternatively, if we superpose the two waves with a pi radian phase difference, then the outcome is a wave of amplitude 
zero with zero energy! The solution lies in the importance of consideration of phases. 
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Figure 4. Flowchart of shape optimization process. We start by specifying the working conditions 
of the device e.g. operating depth, sea bottom depth and incident spectrum, as well as Genetic Al- 
gorithm (GA) parameters. Next, by using our shape generation method, shapes are generated and 
meshed. Then through Boundary Element Method (BEM) solver NEMOH, hydrodynamic coeffi- 
cients of each shape are calculated. Next, the optimum power take-off parameters for each shape are 
found using a separate optimization subroutine. At the next step, the objective function value, i.e. 
the power that can be absorbed by each shape, is calculated. Then the evolutionary process of GA 
optimization continues until reaching to the maximum number of generations where the optimiza- 
tion process stops and the shape with maximum absorbed power is chosen as the optimum shape. 
For each optimization case the maximum number of generations is tuned in a way that the increment 
in absorbed power by the fittest individual in the consecutive generations be small and when there is 
no significant improvement in the values of fitness of the population from one generation to the next. 
On average an optimization case for a broadband directional wave climate required 21.71 hours of 
computation on 8 CPU cores. 


to increase the convergence efficiency of genetic algorithm, e.g. how many parents survive to the next generation, or 
putting a limit on the life of each member of the population. But eventually, after enough number of generations (that 
may range from few to billions or higher) in many cases a global optimum is reached. Genetic algorithm is a heuristic 
search that starts from a random search, but employs an evolutionary logic to find the optimum solution. Genetic 
algorithm is typically computationally expensive, but in many cases its performance can significantly exceed that of 
gradient-based methods, particularly when the objective function is not provided in an analytical form. 

Here the optimization objective function is to maximize P r defined according to (2.1), and variables to be optimized 
are a n ,<p n ,p and q given in equations (2.4) and (2.5) under discussed constraints. The optimization procedure is as 
follows: we initialize our Genetic Algorithm with random selections of a n ,(j) n ,p and q. Each set corresponds to a 
unique shape whose hydrodynamic coefficients can be calculated via NEMOH after proper meshing is generated. 
Since a wide range of different geometries are tested, a careful convergence test for each shape is performed through 
refining the mesh and re-calculating hydrodynamic coefficients. We then using Genetic Algorithm optimize power 
take-off parameters c, k for the shape to yield the maximum power. This is done for every individual shape in the 
population. Then objective function is evaluated and subsequently the next generation is created. The procedure 
continues until convergence is reached. 
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3. Results and Discussion 
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Figure 5 . Optimum absorber plate shapes for monochromatic unidirectional incident waves arriv- 
ing from the right (i.e. 0 = 0°). N c = 3(a), 4(b), 5(c), 6(d) and 7(e), for normalized area A/A 2 = 0.2, 
kliop r = 1.13 and kd = 0.11. For comparison of shapes, we show a superposition of all optimum 
shapes in figure f. Shown in this figure are the equivalent circle of the same area (red dashed-line) 
and optimum shapes (blue solid line) and coefficients of shape generation Fourier modes are given 
in Appendix II. 


We first consider the case of monochromatic unidirectional incident wave, for which we already expect that the 
optimum shape is long and narrow (i.e. stretched along y-axis). Therefore, this step also serves as a validation of 
our approach. Since typically the cost of ocean objects are approximated proportional to the material used, in our 
optimizations we assume that the total area of the device A is constant. If the geometry is a simple circle, we call it the 
equivalent circle whose radius is r e = yjA/n. 

In practice, the range of variables to be optimized is always limited by design constraints that include for example 
fabrication limitations and force distribution and concentration. Here we assume 0.5 < p,q <2 which physically 
means that the optimum shape cannot be too narrow or long. Also for the equivalent circle radius r e and base circle 
radius ro (cf. equation 9) we assume ro/r e = 0.95, a n /ro < 0.40 in order to keep the shape near a circle and prevent 
formation of complex geometries with sharp corners. Clearly the methodology is very general, and we have checked 
that the qualitative nature of results presented here do not change by varying these constraints within a reasonable 
range. 

For a total area of A/A 2 = 0.2, where A is the wavelength of incident unidirectional monochromatic waves, we 
compare the optimum shape (shown in blue solid lines) that gives the maximum normalized average power P r , along 
with the equivalent circle (shown in red dashed-line) in figures 5a-e. We consider, respectively, the number of modes 
N c = 3-7 in figures 5a-e, and an overlap of all shapes in figure 5f. Clearly the optimum absorber plate shape tends 
toward an elongated shape perpendicular to the direction of incident wave (i.e. 0=0), and is symmetric about that axis. 
With N c =3 we obtain a power production P,= 3.29 with the optimum shape depicted in fig. 5a, that is, the optimum 
shape has more than 300% higher energy capture capability than a circular absorber of the same area. By increasing 
the number of modes N c , at the cost of more curves and corners the performance consistently increases to P r = 3.41, 


Table 1. The shape difference factor defined according to (3.1), for comparison of relative 
absorber plate shapes presented in figure 5. Differences in shape decrease as N c increases that 
indicates a convergence of our scheme with N c . 

<j?3,4 <^4.5 <^5,6 <^6,7 

0.47 % 0.46% 0.39 % 0.32% 

<^3,7 <jj4,7 <j?5.7 *?6.7 

3.72% 2.44% 1.77% 0.32% 


3.48, 3.51 and 3.52 for respectively N c = 4, 5, 6 and 7. Clearly the performance is asymptotically approaching an 
upper limit of P r ~ 3.5 with the increase of N c . 

We can also show that shapes converge as N c increases. To do so, we define a shape difference factor £,a,b, as a 
metric to measure the difference between two absorber plate shapes A, B, according to 


(3.1) 


^ ,l = 2n 


r2n 

r B {0)-r A (d) 

Je = o 

r A(0) 


dO. 


We show the values of | for (i) successive shapes, and (ii) each shape compared to the one with highest N c (= 7) in 
table 1. The difference between shapes decreases from 0.47% for N c = 3, 4 to 0.32% for N c = 6, 7, and also the 
difference with each shape and N c = 7 shape also decreases mono tonic ally. Therefore, the optimization is convergent 
with respect to N c , though as it is seen the rate of convergence is relatively slow. 

It can be shown that optimum shapes in figure 5 under mentioned constraints have minimum or very low second 
moment of area about their y-axis (i.e. an axis parallel to the v-axis that goes through each shape’s center of mass). 
This is equivalent of having the narrowest shape within the pool of shapes under above assumptions, which is expected 
since for monochromatic unidirectional waves under potential flow theory in fact the narrowest shape has the highest 
capture width ratio. We would like to comment that in some cases, for instance in the case of figure 5-c, the optimum 
shape has a slightly (< 3%) higher second moment of area than the shape with lowest moment of inertia with N c = 5. 
The reason is that under considered constraints and for a given N c (for the case of figure 5-c, N c = 5), the optimum 
shape for P, as well as the shape with the minimum second moment of area have (smooth) dents and bulges (see 
e.g. figure 5-e). Once these dents and bulges are introduced, second moment of area may not suffice to describe the 
optimum shape and a full wave-analysis is needed to investigate wave-body interaction. That’s the reason why in the 
presence of dents and bulges that appear for higher N c ’s, optimum P, shape may be slightly different than the shape 
with the minimum second moment of area. 


Monochromatic directional incident wave 

Here we consider the case of proposed absorber working in a monochromatic directional sea state, i.e. when all 
waves have the same frequency, but come from different directions. As discussed before (see also Appendix I), in this 
case the relative phases of waves arriving from different direction play an important role. Let’s first consider the case 
in which all waves are in phase with respect to the origin of the coordinate system. We show in figure 6 the optimum 
absorber shape for P r , as calculated from (2.2), for zero directionality of incident waves (figure 6a, which is the same 
as the case in figure 5c), and cases of a directional incident wave with energy evenly distributed across |0| < 20° 
(figure 6b), 0| < 40° (figure 6c), |0| < 60° (figure 6d), and |0| < 80° (figure 6e). Based on convergence behavior for 
the case of unidirectional waves (cf. figure 5), in the following simulations we choose N c = 5. 

In the absence of directionality (figure 6a) as discussed before, the optimum absorber shape has a narrow geometry. 
As the spreading angle departs from zero, say for |0| < 20° in figure 6b, the optimum shape develops two tilted 
wings in order to capture energy from all major incident directions. In other words, optimum shape takes a symmetric 
inclined form in order to increase the length of the plate perpendicular to the incoming waves. This trend is further 
highlighted in the case of 1 0 < 40° (figure 6c). 

For a higher directionality angle of 1 9 \ < 60° (figure 6d), the absorber develops two wings on each side (total of four 
wings). It appears that development of these four wings is to capture energy from even more directions: each wing 
captures the incoming energy normal to its axis, hence adds to the absorber’s capturing orientational diversity. As the 
directionality angle further increases (figure 6e), the optimum absorber shape deviates further from the slender shape 
and tends more toward a circle. This is in fact expected since in the limit of |0| < 180°, waves coming symmetrically 
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Figure 6. Optimum absorber plate’s shapes for monochromatic directional incident waves whose 
energy is uniformly distributed across the provided range of 6 (cf. figure3.a). As the spreading angle 
increases from 0 = 0°, the optimum shape leans symmetrically backward to face maximum area of 
itself normal to the incident wave energy. With a further increase in the directorial range, i.e. for 
1 0 1 < 60°, 80°, a double- wing topology is obtained, as well as a bulge-shaped area in the middle that 
captures energy of high-angle incident waves. In the limit of 1 0 1 < 180° (figure 6f), the optimum 
shape is a circle. Parameters and line conversion is the same as in figure 5 and coefficients of shape 
generation Fourier modes are given in Appendix II. 


from all directions, the optimum shape becomes a circle, as is seen in figure 6f. The trend toward a circle is also seen 
from the development of a bulge near the center that grows with the increase in the spreading angle (see figures 6d, e). 

The normalized absorbed power, as expected, decreases as the spreading angle increases. This is also expected 
since for the case of unidirectional incident waves crests are aligned and optimum shape is simply the narrowest 
shape, whereas for directional spreading case this is not the case and therefore overall efficiency is negatively affected. 
The overall power absorption decreases from P r = 3.48, to respectively P r = 2.58, 2.04, 1.43 and 1.16 as spreading 
angle increases from 0 to 80 degrees (figures 6). It is worth nothing that for incident waves with directional symmetry 
investigated here, as expected all optimum shapes are also symmetric with respect to the mean wave direction (i.e. 

e m = 0 °). 

If the monochromatic directional spectrum arrive with random phases, then the optimum shape depends on those 
specific phases. That is, for every set of random phases, a different optimum shape is obtained. Let’s consider the 
case of figure 6d, in which waves are arriving with directional spread of 1 0 1 < 60°. The optimum shape for eight set 
of random phases (with uniform distribution) is shown in figures 7a-h. The absorbed power in all cases are higher 
than a circular-shape absorber from 10% to 38%, with an average of 25% increase in the efficiency. The optimum 
shape for an unknown set of random phases (which is the case in real ocean), therefore, will be the average of the 
obtained geometries (figure 7i). Calculation of standard error shows a relatively good convergence to the mean (figure 
7j), though shapes for different random phases may deviate from this mean by an average standard deviation of less 
than five percent (see figure 7k). We finally show in figure 71 superposition of all eight geometries and the average 
shape. It is certainly of interest to see how the average geometry (figure 7i) performs under different random phases 
of figure 7a-h. Our simulations show that if the average profile is chosen, under the eight wave conditions discussed 
above, we obtain respectively P, = 1.09, 1.08, 1.17, 1.19, 1.12, 1.25, 1.19 and 1.02, which is on average 14% increase 
in the efficiency compared to a circular-shape geometry. 
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Figure 7. Shape optimization of the absorber plate under monochromatic directional incident 
waves (|0| < 60°) with random phases (c.f. figure 6d which is the same optimization but with 
zero phases). Figures a-h show the optimum shape and normalized power P r for eight different 
sets of random phases (uniformly distributed random numbers from [0, 2zr] ). Clearly, for each set 
of random phases a different optimized shape is obtained. The average enhancement in the power 
extraction is ~ 25%. The average of these shapes is shown in figure (i), along with standard error 
(figure j) and standard deviation (figure k). The average enhancement in the power extraction using 
the average shape is ~ 14%. Finally, we show the superposition of all shapes and the average in 
figure (1). Parameters and variables are the same as in figure 6 and coefficients of shape generation 
Fourier modes are given in Appendix II. 


Similar results are obtained for other angular distributions. For instance, for uniformly distributed monochromatic 
waves across 1 8 | < 40° (cf. figure6c) but with random phases, case-by-case optimization yields an average increase 
in power of 60%, and the average profile obtains an average increase of 42%. Note that the increase in the average 
absorbed power is expected as the angular distribution decreases due to a similar mechanism discussed before. 

Polychromatic (broadband) unidirectional incident wave 

If the incident wave is polychromatic and unidirectional, similar to the case of a monochromatic unidirectional 
wave, the optimum shape must be the most elongated shape. Our optimization obtains identical shapes to those 
reported in figure 5, for each N c . The absorbed power, however, as expected will be much less. This is due to the fact 
that our linear power take-off parameters can only be tuned to one frequency, and therefore the rest of the spectrum 
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is detuned to the power take-off (It should be noted that the optimum power take-off parameters are found through 
optimization as well to ensure maximum extracted power for each shape). For an incident JONSWAP spectrum of sea 
state five (rough sea conditions) with H s = 3.25 meters and T p = 9.7 seconds and for N c 3 7 we obtain, respectively, 
P r = 1.41, 1.52, 1.56, 1.58 and 1.61. The pattern is similar to the case of figure 5, but it is converging to P r ~1.6, i.e. 
our optimum shape performs about 60% better than a circular absorber. 


Polychromatic (broadband) directional incident wave 

For the case of polychromatic directional incident waves, similar to monochromatic directional waves, phases of 
waves play an important role both in the optimum shape and the absorbed power. This is due to the fact that when 
phases are randomly distributed, it is possible that some waves at specific directions undergo a destructive interference 
with other waves of the same frequency coming from another direction, and therefore their effective contributions to 
power absorption decreases. Alternatively, some may undergo constructive interference whereby contribute positively 
to the power absorption. 



FIGURE 8 . Optimization of the absorber shape under the action of polychromatic (broadband) di- 
rectional incident waves, (a-h) optimum shape with eight different sets of random phases, (i) the 
averaged profile, (j) standard error of averaging, (k) standard deviation of the shape distribution and 
(1) optimum shape under the assumption of the entire spectrum being in phase (all phases set to 
zero). Parameters are k p h opr = 1.13, k p d = 0.11, 0.2 [Hz] < / < 4 [Hz], 0 < 80° and coefficients 
of shape generation Fourier modes are given in Appendix II. 

To demonstrate this case, we consider an incident JONSWAP spectrum of sea-state five given by equations (1.4) 
and (1.5) with H s = 3.25 meters and T p = 9.7 seconds, and spreading parameters of ii = 0, ,y = 2, 0 = 160°. For 
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an absorber of area A/X p = 0.064, figure 8a-h shows optimum shapes for eight different sets of random phases. The 
obtained optimum shapes, similar to the case of monochromatic directional waves in figure 7, are different from one 
case to the other, with the two wings more highlighted in almost all figures when compared to the cases shown in 
figure 7a-h. The enhancement in power absorption ranges from 208% to 637%, that is two to more than six times 
higher than a circular absorber of the same area. The average absorption enhancement is 371%. 

We also show the average profile of the absorber (figure 8i), standard error which is less than 1% (figure 8j), and 
standard deviation (figure 8k). The average absorber geometry of figure8i, under each of incident waves of figure8a-h 
obtains respectively P r = 4.94, 1.91, 5.21, 1.86, 3.87, 1.80, 1.94 and 3.89, whose average is P,= 3.17. This means 
that the average profile shown in figure8i has more than 300% higher absorption capability than a circular shape of 
the same area. For the sake of comparison, we also calculate the optimum shape under the assumption of all waves 
being in phase (i.e. we set all phases equal to zero). Figure 81 shows the optimum shape obtained (blue solid line) and 
compares it with the average shape (black dashed line) and the same-area circle (red dashed line). The relative power 
in this case has the surprisingly high value of P r = 7.53, and it is interesting to note that the obtained shape under 
in-phase assumption is very close to the mean shape. 

As discussed before, the methodology developed here does not assume the incident waves having angular symmetry 
with respect to 9 = 0°. It is, therefore, worth investigating a case in which the incident wave spectrum is not symmetric. 
We consider an asymmetric directional JONSWAP spectrum of sea state five (rough sea conditions) with H s = 3.25 
meters and T p = 9.7 seconds, and spreading parameters of jj. = 0.5, © = 160°, and s = 15 (figure 3c). The procedure 
is similar as before, and we consider in-phase incident waves (as in figure81). The optimum shape is shown in figure 
9 for which P r = 4.55. Clearly the shape is asymmetric as a response to an asymmetric incident spectrum. Note that 
the incident spectrum (figure 3c) has an inclination angle of 0 ~ 10° ~ 25°, and accordingly the optimum shape has 
almost the same inclination angle in its vertical axis in order to position itself nearly perpendicular to the direction of 
waves in order to maximize the absorbed power. 



Figure 9. Optimum geometry of the absorber plate under the action of a polychromatic directional 
incident wave whose directionality is not symmetric (cf. figure 3c). Chosen parameters are k p h opr = 
1.13, k p d = 0.11, 0.2 [Hz] < / < 4 [Hz], |0| < 80° and all waves are assumed to be in-phase (cf. 
figure 81) and coefficients of shape generation Fourier modes are given in Appendix II. 


4. Conclusion 

Harnessing energy of water waves is based on conversion of the energy within the waves to drive a power take- 
off unit by means of an intermediate absorber or interface. However, for having efficient energy conversion, the 
wave energy converter must be optimized in the design stage. We presented here a robust and systematic method for 
optimizing the absorber shape of a submerged planar pressure differential wave energy converter, using the Genetic 
Algorithm aiming to improve its wave power absorption level. A new parametric description of absorber shape based 
on Fourier decomposition of geometrical shapes was introduced. For each shape configuration we optimized the shape 
parameters and power take-off characteristics. The shape optimizations were done for different sea states as well 
as a general directional wave spectrum with asymmetric angular distribution. An analytic expression for calculation 
of absorbed power by wave energy converters subject to directional spectrum of waves was presented. Optimum 
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shapes were found to be elongated perpendicular to the mean direction of incident waves with round corners. Having 
symmetry with respect to the mean wave direction for waves with symmetric angular distribution, the optimum shapes 
became asymmetric for an asymmetric angular distribution of incident waves. We showed that the optimum shape 
may have a significant higher energy capturing capability, sometimes nearly an order of magnitude, when compared 
to a circular absorber shape of the same area. 

The focus of the current manuscript is on a device with one degree of freedom (i.e. heave motion). A natural 
extension would be to relax other degrees of freedom and find the optimum shape and corresponding power capture. 
While this is straight forward and follows closely the optimization methodology presented here, it is computationally 
a lot more expensive than a single degree of freedom optimization. It is to be noted that while relaxing more degrees 
of freedom may offer advantages in engineering and fabrication of the wave energy device, it may not necessarily 
result in enhancement of the efficiency. Current manuscript also considers small amplitude incident waves for which 
linear theory is valid. At such a limit, and for scales at which wave energy conversion devices work in high Reynolds 
numbers, viscous effects are known to be of a negligible significance [67]. However, when incident wave amplitude 
is large (e.g. during a storm), nonlinear effects and viscous effects both become important and must be taken into 
account. Specifically, if waves are large enough that they break over the device, then certainly both nonlinearity and 
viscosity play major roles in the dynamics of wave-device interaction. In this case, direct simulation of full governing 
equation is inevitable. 


Appendix I: Energy damping and the superposition principle 

The objective here is to present expressions to calculate generated (damped) energy in a damper of a linear mass- 
spring-damper system under the action of two sinusoidal forcing functions. In brief, if the two forcing functions 
have different frequencies, then generated energy is the linear superposition of energy generated under each of the 
excitation components. Nevertheless, if the two excitation forces have the same frequency (e.g. two same-frequency 
waves arriving from different directions), then the relative phases of the two waves play a critical role in the overall 
energy production. Overall energy production in this case may become zero in the special case of the two phases being 
7T-Radian different (at the location of the device). The energy production is maximum if the two phases are the same. 
While this is a classic subject, we repeatedly find mistakes and confusion in the work of researchers, and therefore 
decided to discuss this here. 


F, F 2 F,+F 2 



a) b) c ) 


Figure 10. Mass-spring-damper system subject to (a) external force with frequency 
C0 \ ; (b) external force with frequency © 2 ; and (c) superposition of cases (a) and (b) 


Consider a linear mass-spring-damper system under the action of sinusoidal forces F\ ( C0j ) = F\ = 

P 2 e'( a>2t+ ‘I >F ’ 2 \ and the linear superposition of the two (figure lOa-c), where (Oj,(f)Fj are frequency and phase of each 
excitation. Under the action of each of the single forcing, corresponding system response is zj = with 

(4.1a) zj = ■ Fj , 

y (K — Ma>j) 2 + (C C0j) 2 

(4. lb) = <j>Fj - 0k, c = <t>F,j - arctan 

Therefore, the absorbed power is 

(4.2) P j (t)=Re{F c j(t)}Re{Vj(t)}=C{Re{z j (t)}) 2 =C[-(O j Zjsm((Ojt + 0 z j)] 2 , 
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C(Qj \ 
K — Mcoj )' 



where F c j is the force of damper and Vj is the vertical velocity of the mass M due to external force Fj. Therefore, if 
the system of mass-spring-damper is subjected to two forcings, then 

(4.3) P(t) = C [-(Oizi sin ( o>i f + 0- 1 ) - Ohzi sin(ohf + 0,-,2)] 2 - 

The quantity of interest, however, is the time-averaged power 

(4-4) = ^ [ T P(t)dt, 

1 Jo 

where T is integration time, which for single frequency cases is the period of the external forcing, and for multiple fre- 
quency cases is integer multiples of Least common multiple of periods (if does not exist, then it is infinity). Therefore, 
it is always accurate to calculate Pf ve in the limit of T — > °°. 

For cases of interest discussed above, assuming that a>i ^ ah the average power becomes 

11C 

(4.5) Pj lve = -Ccofzi + - + -(bounded value), 
and if ffli = (So 

1 C 

(4.6) Pj ve = -Cm\ (z\ + zl) +Ca>f ziZ 2 Cos(</>i - <j> 2 ) + -(bounded value). 

Clearly the last term in both expressions go to zero as T — > °°. Therefore, effect of relative phases of the two forcing 
excitations with different frequencies does not affect the overall power production, but if they have the same frequency 
then these phases play an important role. Expressions (4.5) and (4.6) (without the last term) are used to calculate total 
energy in different cases studied in this manuscript. 

Appendix II: Optimum coefficients 

In the following, are shape coefficients (cf. equations 9, 10) where a nr = a n /ro , ro is the base circle 

radius, and (0/(0 n are non-dimensional parameters. 

Table 2. Optimum shape coefficients for Fig. (5); £ = C/[2y/K(m+A)\, (o n = K/(m+A ) 
where m,A,C,K are absorber’s mass and added mass, and PTO’s optimal damping and restoring 
coefficients respectively 
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Table 3. Optimum shape coefficients for Fig. (6); C, = C/[2^/K(m+A )\ , a>„ = \J K/(m+A ) 
where m,A,C,K are absorber’s mass and added mass, and PTO’s optimal damping and restoring 
coefficients respectively 
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Table 4. Optimum shape coefficients for Fig. (7); £ = C/[2y/K(m+A)\, co n = K/(in+A ) 
where m,A,C,K are absorber’s mass and added mass, and PTO’s optimal damping and restoring 
coefficients respectively 
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Table 5. Optimum shape coefficients for Fig. (8); £ = C /\2\J K(m)\, a>„ = \J K /(m) where 
m,C,K are absorber’s mass, and PTO’s optimal damping and restoring coefficients respectively, 
c o p is the angular frequency of spectrum’s corresponding peak amplitude 
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Table 6. Optimum shape coefficients for Fig. (9); £ = c/\2\Jk{m)\, a> n = \Jk]Jrn) where m,c,k 
are absorber’s mass, and PTO’s optimal damping and restoring coefficients respectively, co p is the 
angular frequency of spectrum’s corresponding peak amplitude 
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Appendix III: Derivation of expression for the area of the optimization candidate shapes 

For the shapes built using equations (2.4) and (2.5) it was mentioned in section 2 that the area of the shapes can be 
found by equation (2.6) which by using Green’s theorem [65] here we derive it analytically. 

Green’s theorem asserts that if 

(4.7) y: 0 ha (x(0),y(0)) (0< 0 <2n) 
is a closed curve bounding counter-clockwise a region B C R 2 , then 

1 r 2n 

(4.8) Area(B) = - / \x(0)y(0) —y(0)x(0)]d0. 

2 Jo 

Now, for a shape generated using the equations below 

N c 

(4.9) r(0) = r o +£a„sin(n0 + 0„), 

n = I 

(4.10) x = pr(9) cos(0), y = qr(9) sin(0), 
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by replacing the expressions of x.y into equation (4.8) and calculating the integral by hand or a software the area of 
the shape becomes 


(4.11) 



As a check one can put p = q = 1 and do not add Fourier modes to the circle (i.e. a n = 0) in order to recover the known 
formula for area of a circle with radius ro as nrfy 
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